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, ABSTRACT 

in ; 

f*^ > Aims. We present 3.65" x 3.34" angular-resolution IRAM Plateau de Bure Interferometer (PdBI) observations of the CS 7=2—1 line toward 
. the Horsehead Photodissociation Region (PDR), complemented with IRAM-30m single-dish observations of several rotational lines of CS, 
' C 34 S and HCS + . We analyse the CS and HCS + photochemistry, excitation and radiative transfer to obtain their abundances and the physical 
^ (— i conditions prevailing in the cloud edge. Since the CS abundance scales to that of sulfur, we determine the gas phase sulfur abundance in the 
Q | PDR, an interesting intermediate medium between translucent clouds (where sulfur remains in the gas phase) and dark clouds (where large 
depletions have been invoked). 

Methods. A nonlocal non-LTE radiative transfer code including dust and cosmic background illumination adapted to the Horsehead geometry 
has been developed to carefuly analyse the CS, C 34 S, HCS + and C 18 rotational line emission. We use this model to consistently link the line 
observations with photochemical models to determine the CS/HCS + /S/S + structure of the PDR. 

Results. Densities of n{H2) — (0.5 - 1.0)xl0 5 cirT 3 are required to reproduce the CS and C 34 S 7=2-1 and 3-2 line emission. CS 7=5-4 lines 
show narrower line widths than the CS low-7 lines and require higher density gas components not resolved by the ~10" IRAM-30m beam. 
These values are larger than previous estimates based in CO observations. We found ,^(CS )=(7±3)xl0~ 9 and / ^(/7CS + )=(4±2)xlO~" as the 
averaged abundances in the PDR. According to photochemical models, the gas phase sulfur abundance required to reproduce these values is 
d ' S/H=(3.5±1.5)xl0~ 6 , only a factor <4 less abundant than the solar sulfur elemental abundance. Since only lower limits to the gas temperature 
are constrained, even lower sulfur depletion values are possible if the gas is significantly warmer. 

Conclusions. The combination of CS, C 34 S and HCS + observations together with the inclusion of the most recent CS collisional and chemical 
rates in our models implies that sulfur depletion invoked to account for CS and HCS + abundances is much smaller than in previous studies. 

Key words. Astrochemistry - ISM clouds - molecules - individual object (Horsehead nebula) - radiative transfer - radio lines: ISM 
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1. Introduction 

Sulfur is an abundant element (the solar photosphere abun- 
dance is S/H=1.38xl0~ 5 ; |Asplund et al. 2005|l, which remains 
undepleted in diffuse interstellar gas (e.g. iHowk et al. 2 006) 
and HII regions (e.g. Martin-Herna ndez et al. 2 002 
|Garcia-Rojas et al. 2006| and references therein) but it is 
historically assumed to deplete on grains in higher density 
molecular clouds by factors as large as ~10 3 (Tieftrunk et al. 
1994). This conclusion is simply reached by adding up the 
observed gas phase abundances of S-bearing molecules in 
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well known dark clouds such as TMC1 (e.g. Irvine et al. 1985; 
Ohishi & Kaifu 1998). As sulfur is easily ionized (ionization 
potential ~ 10.36 eV), sulfur ions are probably the dominant 
gas-phase sulfur species in translucent gas. Ruffle et al. (1999) 
proposed that if dust grains are typically negatively charged, 
S + may freeze-out onto dust grains during cloud collapse 
more efficiently than neutral species such as oxygen. However, 
the nature of sulphur on dust grains (either in mantles or 
cores) is not obvious. Van der Tak et al. (2003) observed 
large abundances of gas phase OCS, ~10~ 8 , in star forming 
regions, and suggested that together with the detection of solid 
OCS (with an abundance of ~10~ 7 ; Palumbo et al. 1997), 
it implies that OCS is a major sulfur carrier in dust grains. 
However, the ~4.9 fim ice feature attributed to OCS is best 
reproduced when OCS is mixed with methanol. In fact, the 
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band is blended with a methanol overtone whose contribution 
has not been studied in detail (Dartois 2005). In any case, the 
absence of strong IR features due to S-bearing ices in many 
ISO's mid-IR spectra (e.g. Boogert et al. 2000; Gibb et al. 

2004) and the presence of S n recombination lines in dark 
clouds such as Rho Ophiuchi (Pankonin & Walmsley 1978) all 
argue against a large depletion of sulfur from the gas phase. In 
this case, the abundance of species such as CS may indicate 
that something important is lacking from chemical models 
or that an abundant sulfur-bearing carrier has been missed. 
Therefore, the abundances of sulfur species remain interesting 
puzzles for interstellar chemistry. In the case of dense clouds, 
standard chemical models predict that most of the gas phase 
sulfur is shared between S, SO and CS (Millar & Herbst 1990), 
while H2S is also abundant in the Orion Bar PDR (Jansen et 
al. 1995). In all these cases, a large sulfur depletion, ~10 2 , was 
required in the models to explain the observed abundances. 

PDRs offer an ideal intermediate medium between dif- 
fuse and dark cloud gas to investigate the sulfur depletion 
problem. In this work we have tried to determine the CS 
abundance in the Horsehead PDR as a tool for estimat- 
ing the sulfur gas phase abundance. However, CS chem- 
istry is an open issue itself in different environments, from 
hot cores (e.g. Wakela m et al. 200 4 1 to extragalactic sources 
(e.g. Ma rtin et al. 20 05 1. Recent laboratory experiments on dis- 
sociative recombination of HCS + and OCS + (Montaigne et al. 

2005) imply a substantial modification of previous reaction rate 
coefficients, dissociative channels and branching ratios used in 
chemical models. The latest available reaction rates and col- 
lisional coefficients have been used in our photochemical and 
radiative transfer models. 

1.1. The Horsehead nebula 

The Horsehead nebula, appears as a dark patch of ~5' di- 
ameter against the bright HII region IC434. Emission from 
gas and dust associated with this globule has been detected 
from IR to millimeter wavelengths fAbergel et al. 2002 
|Abergel et al."2003l IPound et al. 2"003l |Teyssier et al"2004 
Habart et al. 2005, Pety et al. 2005a), although the first 
astronomical plates were taken ~120 yr ago. In particular, the 
Horsehead western edge is a PDR viewed nearly edge-on and 
illuminated by the 09.5V star crOri at a projected distance 
of ~3.5pc ( |Anthony-Twarog 1982) . The intensity of the 
incident FUV radiation field is% ^60 relative to the interstellar 
radiation field (ISRF) in Draine's units dDraine 1 978 1. 

According to the evolutionary view of Reipurth & Bouchet 
(1984), the Horsehead nebula was a quiescent and dense cloud 
core embedded in a more diffuse cloud (L1630). The erosive 
action of the UV radiation from crOri on the ambient gas led 
to the apparent emergence of the core cloud, as in the earli- 
est stages of Bok globules still attached to their parental cloud. 
However, the observed morphology together with the veloc- 
ity gradients of the cloud, require a more involved description 
including a pre-existing rotating velocity field as well as den- 
sity inhomogeneities in the initial structures (Pound et al. 2003; 
Hily-Blant et al. 2005). The erosive effect of the ionizing and 



Table 2. Line parameters for the IRAM 30-m CS observations. 



Molecule 


Transition 


Frequency 


HPBW 






(GHz) 


(arc sec) 


CS 


7=2-1 


97.980968 


25 




7=3-2 


146.96905 


16 




7=54 


244.93561 


10 


c 34 s 


7=2-1 


96.412982 


25 




7=3-2 


144.61715 


16 


HCS + 


7=2-1 


85.347884 


29 



dissociating radiation field together with these initial condi- 
tions explain the peculiar shaping of the Horsehead nebula. In 
particular, the densest regions of the initial inhomogeneities are 
now believed to be the East- West filamentary material connect- 
ing it to the parental cloud, and the PDR. In this work we have 
studied the PDR through CS, C 34 S and HCS + observations. 

2. Observations and data reduction 

2.1. Observations 

2.1 .1 . Pico Veleta single-dish 

The single-dish data presented in this paper have been gath- 
ered between February and October 2004 at the IRAM 30-m 
telescope. The Horsehead nebula PDR was mapped in the CS 
7=2-1 and 5—4 lines in order to provide the short-spacings for 
the interferometric observations presented thereafter. The final 
map consists of 5 on-fhe-fly coverages performed along per- 
pendicular scanning directions, and combined with the PLAIT 
algorithm introduced by Emerson & Grave (1988), allowing to 
efficiently reduce the stripes over the map. The noise levels (lcr 
rms) per regridded pixel and resolution channel of 80 kHz are 
of the order of 0.15 K at 3 mm, and 0.64 K at 1.3 mm. The lat- 
ter value was not low enough to provide any useful mapping 
information at 1.3 mm since the CS 7=5-4 line peak is <1 K. 

In complement to these data, dedicated positions were 
probed over a larger set of species and transitions. The fre- 
quency switching mode was used to observe CS 7=2-1, 3-2 
and 5-4 lines, as well as C 34 S 7=2-1, 3-2, and HCS + 7=2- 
1 lines. Table |2 summarizes the corresponding observing pa- 
rameters. Longer integrations allowed to reach, in a resolution 
channel of 40 kHz, rms noise levels of 25, 42 and 36 mK at 3, 2, 
and 1.3 mm respectively. All CS and C 34 S lines were detected 
with a S/N ratio better than 10. Figs.|4]and[8]show some spectra 
collected at positions inside and across the PDR. 

The data were first calibrated to the scale using the so- 
called chopper wheel method (Penzias & Burrus 1973), and 
finally converted to main beam temperatures using efficiencies 
(Beff/F eff ) of 0.8 1 , 0.74 and 0.50 at 3, 2 and 1 .3 mm respectively. 

2.1.2. Plateau de Bure Interferometer 

PdBI observations dedicated to this project were carried out 
with 6 antennae in the BCD configuration (baseline lengths 
from 24 to 331 m) from August 2004 to March 2005. The 
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Fig. 1. Integrated emission maps obtained with the Plateau de Bure Interferometer. The center of all maps has been set to the 
mosaic phase center: RA(2000) = 05h40m54.27s, Dec(2000) = -02°28'00". The map size is 110" x 110", with ticks drawn 
every 10". The synthesized beam is plotted in the bottom left corner. The emission of all lines is integrated between 10.1 and 
11.1 km s" 1 . Values of contour level are shown on each image wedge. The four panels are shown in a coordinate system adapted 
to the source: i.e. maps have been rotated by 14° counter-clockwise around the image center to bring the exciting star direction 
in the horizontal direction as this eases the comparison of the PDR models. Maps have also been horizontally shifted by 20" to 
set the horizontal zero at the PDR edge, delineated by the vertical line. 



580 MHz instantaneous IF-bandwidth allowed us to simulta- 
neously observe CS, I-C3H and 34 SO at 3 mm using 3 different 
20MHz-wide correlator windows. Another window was cen- 
tered on the 13 CO 7=2-1 line frequency at 1 mm. The full IF 
bandwidth was also covered by continuum windows both at 3.4 
and 1 .4 mm. Only CS 7=2-1 and 13 CO 7=2-1 (not shown here) 
were detected. 

We observed a seven-field mosaic. The mosaic was 
Nyquist sampled in declination at 3.4 mm and Nyquist sampled 
in Right Ascension at 1.3 mm. This ensures correct sampling in 
the illuminating star direction both at 3 and 1 mm while max- 
imizing the field of view along the edge of the PDR eases the 
deconvolution. This mosaic, centered on the IR peak (Abergel 



et al. 2003), was observed for about 30 hours of telescope time 
with 6 antennas. This leads to an on-source integration time 
of useful data of 10 hours after filtering out data with track- 
ing errors larger than 1" and with phase noise worse than 40° 
at 3.4 mm. The rms phase noises were between 15 and 40° 
at 3.4 mm, which introduced position errors < 0.5". Typical 
3.4 mm raw resolution was 2.5". 
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Table 1. Observation parameters. 



Phase center Number of fields 

Mosaic 1 a-2000 = 05''40'"54.27 t d 20 oo = -02°28'00" 7 
Mosaic 2 a< 20 oo = 05''40'"53.00 I d 200 o = -02°28'00" 4 



Molecule & Line 


Frequency 


Beam 


PA 


Noise" 


Obs. date 




(GHz) 


(arcsec) 


C) 


(Kkms- 1 ) 




Mosaic 1 












I2 CS 7=2-1 


97.981 


3.65 x 3.34 


48 


1.2x10-' 


Aug. & Oct. 2004 and Mar. 2005 


C l8 7=2-1 


219.560 


6.54x4.31 


65 


9.8xl0~ 2 


Mar. 2003 


Mosaic 2 












12 CO 7=1-0 


115.271 


5.95 x 5.00 


65 


1.2x10-' 


Nov. 1999 


12 CO 7=2-1 


230.538 


2.97 x 2.47 


66 


1.7x10"' 


Nov. 1999 



8 The noise values quoted here are the noises at the mosaic center (Mosaic noise is inhomogeneous due to primary beam correction; it steeply 
increases at the mosaic edges). Those noise values have been computed in 1 kms~' velocity bin. 



Table 3. Calibrator fluxes in Jy. 



Date 


B0420- 


-014 


B0607- 


-157 




3 mm 


1 mm 


3 mm 


1 mm 


20.08.2004 


3.4 


2.9 


1.4 


0.93 


04.10.2004 


3.4 


2.9 


1.6 


0.90 


27.02.2005 


3.5 


2.9 


1.6 


0.90 


02.03.2005 


3.2 


2.3 


1.6 


0.89 


12.03.2005 


3.2 


2.3 


1.6 


0.90 


13.03.2005 


3.2 


2.3 


1.6 


1.00 



2.2. Data processing 

All data reduction was done with the GILDAS 1 softwares sup- 
ported at IRAM (Pety 2005b). Standard calibration methods 
using nearby calibrators were applied to all the PdBI data. The 
calibrator fluxes used for the absolute flux calibration are sum- 
marized in Table |3 

Following IGueth et al. 19961 single-dish, fully sampled 
maps obtained with the IRAM-30m telescope (see sec- 
tion 2.1.1) were used to produce the short-spacing visibilities 
filtered out by each mm-interferometer (e.g. spatial frequencies 
between and 15 m for PdBI). Those pseudo-visibilities were 
merged with the observed, interferometric ones. For imaging, 
we followed the method described by Pety et al. (2005) to pro- 
duce images in different lines of the same source. This results 
in inhomogeneous noise. In particular, the noise strongly in- 
creases near the edges of the field of view. To limit this effect, 
the mosaic field-of-view is truncated. 

Moreover, the natural synthesized beam (3.58" x 1.89" 
at PA 37°) is very asymmetric due to the low declination 
of the Horsehead nebula. We chooe to taper the weights of 
the uv visibilities before imaging using a Gaussian of size 
400 mx 1 15 m at PA 80° to obtain an almost circular clean beam 
(3.65" x 3.34" at PA 48°). This considerably eased the decon- 
volution as 1) the larger beam increases the brightness sensitiv- 
ity and 2) the secondary sidelobes of the dirty beam are much 
less patchy. Deconvolved image nevertheless still shows low- 
level, fake, patchy structures at the scale of the clean beam, 

1 See http://www.iram.fr/IRAMFR/GILDAS for more informa- 
tion about the GILDAS software. 



Table 4. IRAM-30m line observation parameters from 
Gaussian fits. 



Molecule/ 


Offset 


AVp WHM 


jr A dv 


Transition 


(") 


(km s" 1 ) 


(Kkms-') 


CS 7=2-1 


-52,-40 


0.75+0.01 


2.60+0.01 




-64,+30 


0.89+0.01 


3.63+0.01 




-35,-25 


0.78+0.01 


3.68+0.01 




-20,-15 


0.77+0.01 


3.55+0.01 


CS 7=3-2 


-52,-40 


0.72+0.01 


1.82+0.02 




-64,+30 


0.93+0.01 


2.58+0.02 




-35,-25 


0.76+0.01 


2.73+0.02 




-20,-15 


0.80+0.01 


2.40+0.02 


CS 7=5^1 


-52,-40 


0.43+0.02 


0.35+0.01 




-64,+30 


0.76+0.02 


0.62+0.02 




-35,-25 


0.58+0.02 


0.52+0.01 




-20,-15 


0.60+0.03 


0.40+0.02 


C :,4 S 7=2-1 


-52,-40 


0.47+0.02 


0.26+0.01 




-64,+30 


0.67+0.04 


0.38+0.02 




-35,-25 


0.59+0.02 


0.40+0.01 




-20,-15 


0.58+0.03 


0.45+0.02 


C 34 S 7=3-2 


-52,-40 


0.48+0.04 


0.20+0.01 




-64, +30 


0.74+0.05 


0.28+0.02 


HCS + 7=2-1 


-52,-40 


0.8+0.3 


0.07+0.03 




-35,-25 


0.6+0.3 


0.05+0.03 




-20,-15 


0.9+0.4 


0.10+0.03 



mainly along the vertical direction. This is a known artifact of 
the Hogbom CLEAN algorithm when a large spatial dynamic 
(field-of-view/resolution ~ 110/3.5 = 30) is combined with 
high enough signal-to-noise ratio. 

3. Results 

The PdBI integrated emission map of the CS 7=2-1 line, com- 
plemented with previous maps of CO 7=1-0, 2-1 and C ls O 
7=2-1 lines (Pety et al. 2005a), is shown in Fig.^ Integrated 
emission profiles for these lines and for the small hydrocarbon 
C-C3H2 2i2-loi line emission and 1.2 mm continuum emission 
(Pety et al. 2005a) are shown in Fig. [2 In Fig.^the four pan- 
els are shown in a coordinate system adapted to the source: i.e. 
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Fig. 3. Spectra along the direction of the exciting star at 6y = 0". 12 CO 7=1-0, C ls O 7=2-1 and 12 CS 7=2-1 spectra cubes 
were smoothed by a 15"-FWHM ID-Gaussian along the y direction perpendicular to the illuminating star direction. Due to their 
small field of view (in particular in the y direction), the 12 CO 7=2-1 data were just smoothed by a 5"-FWHM circular Gaussian. 



maps have been rotated by 14° counter-clockwise around the 
image center to bring the exciting star direction in the horizon- 
tal direction as this eases the comparison with edge-on PDR 
models. Maps have also been horizontally shifted by 20" to the 
east in order to set the horizontal zero at the PDR edge (delin- 
eated by the vertical line). Therefore, the coordinate system is 
converted from 6RA to 6x and from 6DEC to 6y in arcsec (see 
Fig [D- In the following, our spatial scale to interpret the PdBI 
maps will refer to the 6x scale. 

PdBI observations show that the CS emission is strikingly 
different from that of other species observed at comparable 
spatial resolutions. CS does not follow the filamentary pattern 
along the cloud edge revealed by mid-IR (Abergel et al. 2003), 
H 2 (Habart et al. 2005) or CCH and c-C 3 H 2 emission (Pety et 
al. 2005a). Instead, the behavior of the CS 7=2-1 line emission 
is more similar to that of CO and shows a steep rise across the 
PDR and a plateau in the shielded region. Besides, the CS 7=2- 
1 line emission peak occurs in the well-shielded regions and 



does not coincide with the C ls O 7=2-1 nor with the 1.2 mm 
continuum peaks (i.e. the temperature weighted density peaks, 
see Fig. |2j. Therefore, PdBI observations suggests that CS is 
more abundant in the lower density regions and/or it shows line 
opacity effects in the denser regions. These results confirm that 
there are differences in the spatial distribution of small hydro- 
carbons (Pety et al. 2005a) and other species with similar exci- 
tation requirements. PdBI CO, C I8 and CS line spectra along 
the direction of the exciting star (Sy = 0") are shown in Fig. [3] 
The CS 7=2-1 emission peak and line widths are comparable 
to those of C ls O 7=2-1. Nevertheless, the CS emission distri- 
bution profile in the 6x direction is less smooth, which can be 
an artifact of the data reduction and/or related to the CS photo- 
chemistry. The CS emission increases like the C 18 emission 
in the PDR edge but in the shielded gas, CS rises when C ls O 
decreases. On the other hand, 12 CO peaks closer to the PDR 
edge and lines have broader line widths (Pety et al. 2005a) in- 
dicative of their much larger opacities. 
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CS and C S single-dish observations at larger spatial 
scales are presented in Fig.0] CS line ratios are similar in all 
observed PDR positions. However, there is a trend for CS lines 
to peak where C 18 emission decreases. Again, this is indica- 
tive of larger abundances in the lower density regions and/or 
line opacity effects in the denser regions. The latter hypothe- 
sis is playing a role because CS 7=3-2 lines show asymmetri- 
cal profiles in the whole region, especially red-wing like self- 
absorptions. See for example the (-52, -40) position with re- 
spect to the IRAM-30m C I8 7=2-1 map of Hily-Blant et al. 
(2005; Fig.gJ. In addition, CS 7=3-2 and 2-1 lines must be 
optically thick since their intensity is only a factor ~5 stronger 
than the analogous C 34 S lines, significantly lower than the 
32 S/ 34 S=23 solar isotopic ratio (Bogey et al. 1981). In addi- 



-7.5" < 6y < 7.5" 




<5x (arcsec) 



Fig. 2. Emission profiles along the exciting star direction (PA 
= -104° in the equatorial coordinate system). To improve the 
signal-to-noise ratio, emission profiles have been integrated 
along the perpendicular direction between -7.5" < 6y < 
+7.5". We show from top to bottom 12 CO 7=2-1, 12 CO 7=1-0, 
C I8 7=2-1, 12 CS 7=2-1, c-C 3 H 2 2i2-l i and 1.2 mm dust 
continuum emission (Pety et al. 2005a). The 3cr noise level is 
indicated by the dashed lines. It rises at the cut edges due to 
the primary beam correction. Note that the fields of view of the 
12 CO and C ls O data are smaller than the field of view of the 
12 CS data because of the smaller mosaic size and/or the higher 
frequency. 



Table 5. Einstein coefficients, transition upper level energies 
and critical densities for the range of temperatures considered 
in this work. 



Molecule 


Transition 






Ei 


n„ 






(s- 1 ) 




(K) 


(cm- 3 ) 


C ls O 


7=2-1 


6.01x10 


-7 


15.8 


~8xl0 J 


HCS+ 


7=2-1 


1.11x10 


-5 


6.1 


~5xl0 4 


c 34 s 


7=2-1 


1.60x10 


-5 


6.9 


~4xl0 5 




7=3-2 


5.79x10 


-5 


13.9 


-1x10 s 


CS 


7=2-1 


1.68x10 


-5 


7.1 


~4xl0 5 




7=3-2 


6.07x10' 


-5 


14.1 


-lxlO 6 




7=5-4 


2.98x10" 


-4 


35.3 


~5xl0 6 



tion, CS line ratios are ~0.7 while the optically thin C 34 S 
l^r line ratios are larger ~0.9. Therefore, the CS 7=3-2 line is 
likely to be the most opaque CS line. Finally, Fig.|S]shows clear 
detections of the HCS + 7=2-1 line. As the expected HCS + 
abundance is lower than that of C 34 S, these lines are weak and 
should be optically thin. Line intensities are quite similar in all 
observed PDR positions. 

4. Numerical methodology 

4.1. Photochemical models 

We have used the Meudon PDR code (publicly available at 
http://aristote.obspm.fr/MIS/), a photochemical model of a 
unidimensional stationary PDR (Le Bourlot et al. 1993). The 
model has been described in detail elsewhere (Le Petit et al. 
2006). In few words, the PDR code solves the FUV radia- 
tive transfer in an absorbing and diffusing medium of gas and 
dust. This allows the explicit computation of the FUV radia- 
tion field (continuum+lines) and therefore, the explicit integra- 
tion of consistent C and S photoionization rates together with 
H2, CO, 13 CO, and C 18 photodissociation rates. Penetration 
of FUV radiation into the cloud strongly depends on dust prop- 
erties through dust absorption and scattering of FUV photons. 
Properties of dust grains are those described in Pety et al. 
(2005). We have taken a single dust albedo coefficient of 0.42 
and an scattering asymmetry parameter of 0.6. 

Once the FUV field has been determined, the steady- 
state chemical abundances are computed for a given 
chemical network. The Ohio State University (osu) gas- 
phase chemical network (osu. 2005; September 2005 release; 
http://www.physics.ohio-state.edu/~eric/research.html) has 
been used as our chemical framework. The most important 
changes compared to previous versions are the decrease, by 
a factor of 2, of rate coefficients of photoionization and pho- 
todissociation reactions produced by cosmic-ray-induced H2 
secondary photons, the inclusion of fluorine (F) and its chem- 
istry (see Neufeld et al. 2005) and the update of several reac- 
tion rates. In addition, several changes have been carried out 
by us on the chemical network. In particular, we have intro- 
duced different I8 bearing species into the chemical network 
by assuming similar reaction rates to those involving the major 
isotopologues. Fractionation reactions have been added follow- 
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Fig. 4. IRAM-30mCS 7=2-1, 3-2 and 5^1, and C S 7=2-1 and 3-2 single-dish observations (histograms) at different positions 
of the Horsehead PDR single-dish C ls O 7=2-1 emission centered at a 2 ooo = 05''40 m 58 ! , £ 2 ooo = -02°27'20" (from Hily-Blant 
et al. 2005). 



ing Graedel et al. (1982) and specific photodissociation cross- 
sections for C ls O are explicitly introduced to compute the cor- 
responding photodissociation rate. When available, we have 
also used the photodissociation rates given by van Dishoeck 
(1988), which are explicitly calculated for the Draine interstel- 
lar radiation field (ISRF). Finally, we have further upgraded 
the sulfur network by adding the most recent reaction rates, 
dissociation channels and branching ratios of HCS + and OCS + 
dissociative recombination (Montaigne et al. 2005) and by in- 
cluding the CS photoionization (ionization potential ~11 eV). 
These processes have direct impact on CS chemistry. The re- 
sulting network involves ~450 species and ~5000 reactions. 
Finally, the model computes the thermal structure of the PDR 
by solving the balance between the most important processes 
heating and cooling the gas (see Le Bourlot et al. 1993). Our 
standard conditions for the model of the Horsehead PDR in- 
clude a power-law density profile (Eq. 2) and a FUV radiation 
field enhanced by a factor x — 60 with respect to the Draine 
ISRF (see tableHJl. Different sulfur gas phase abundances, S/H, 
have been investigated. To be consistent with PdBI CO obser- 
vations, thermal balance was solved until the gas temperature 
reached a minimum value of 30 K, then a constant temperature 
was assumed. 

4.2. Radiative transfer models 

We have used a simple nonlocal non-LTE radiative transfer 
code to model our millimeter line observations. The code 
handles spherical and plane-parallel geometries and accounts 
for line trapping, collisional excitation, and radiative excitation 
by absorption of microwave cosmic background and dust con- 
tinuum photons. Arbitrary density, temperature or abundance 
profiles, and complex velocity gradients can be included. A 
more detailed description is given in the Appendix. The choice 
of a nonlocal treatment is needed to analyze optically thick 
lines of abundant, high-dipole moment molecules, such as 
CS, in regions where the gas density is below the critical 
densities of the associated transitions. Table [5] shows the 
critical densities of observed C 18 0, HCS + and CS lines. Under 



Table 6. Horsehead standard conditions and gas phase abun- 
dances. 



Parameter 


Value 


Radiation fields- 


60 (Draine units) 


Cosmic ray ionization rate £ 


5x ltr 1 V 


Density profile n H = n(H) + 2n(H2) 


50 to 2 x 10 5 cnr 3 


Line of sight spatial depth Idepth 


0.05-0.1 pc 


Line of sight inclination angle <p 


0" to 5° 


RefH=n(He)/n H 


1.00 x 10-' 


O/H 


3.02 x 10~ 4 


C/H 


1.38 x 10~ 4 


N/H 


7.95 x 10~ 5 


ls O/H 


6.04 x 10~ 7 


Cl/H 


1.00 x 10- 7 


Si/H 


1.73 x 10~ 8 


Mg/H 


1.00 x 10~ 8 


F/H 


6.68 x 10- 9 


Na/H 


2.30 x 10- 9 


Fe/H 


1.70 x 10~ 9 


P/H 


9.33 x 10- 10 



these conditions, radiative transfer and opacity effects may 
dominate the line profile formation. Our radiative transfer 
analysis has been used to infer abundances and physical 
conditions directly from observations but also to predict line 
spectra from the photochemical model results. The following 
temperature dependent collisional rate coefficients 2 have been 
adopted: 

- for CS: we have used the latest CS + He collisional 
rates from Lique et al. (2006), kindly provided by F. Lique 
prior to publication, scaled by the reduced mass factor 
(Hcs-h 2 I Hcs-He) 1 ^ 2 ■ Most of the models were repeated with 
the older collisional rates of Turner et al. (1992). 

2 Some of them retrieved from BASECOL, a data base for colli- 
sional excitation data at http://www.amdpo.obspm.fr/basecol. 
We considered H 2 , He and H as the collisional partners in all CS, 
C 34 S, HCS + and C 18 excitation models. See Appendix. 
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- for C 34 S: same as CS but using C 34 S spectroscopy to 
compute collisional excitation rates through detailed balance. 

- for C 18 0: CO + H 2 de-excitation rates from Flower (2001) 
but using C 18 spectroscopy to compute collisional excitation 
rates through detailed balance. 

- for HCS + : HCS + + He collisional rates from 
Monteiro (1984), scaled by the reduced mass factor 
(^Hcs + -H 2 /^HCS*-He) 1/2 , have been used. 

5. Modeling and interpretation 

5.1. CS, C 34 S and HCS + single-dish emission 

In order to get a first order approximation of the CS excitation 
and column density, we have assumed that level populations 
are only determined by a Boltzmann distribution at a single 
rotational temperature. If one accepts that lines are optically 
thin, this approach corresponds to the widely used rotational- 
diagram. However, observed CS/C 34 S intensity ratios, and CS 
line profiles (see Fig.|4j clearly show that the low-/ CS lines 
are optically thick towards the Horsehead. Therefore, we have 
included optical depth effects and built a rotational-diagram 
corrected for opacity through: 



N' hin Ej 
In — — + lnC r = In N - InQ - 

gi Trot 



(1) 



where N' hm are the upper level populations determined from 
observations in the optically thin limit (underestimated if lines 
are optically thick), Ej is the upper /-level energy in K, Q is the 
partition function at T rot and C T is the line opacity correction 
factor Tl Lr tJ >1 (Goldsmith & Langer 1999). We have per- 
formed this analysis at different cloud positions. Resulting dia- 
grams are shown in Fig.|5]as a function of different CS 7=2-1 
line opacities (j2-\ — 0, 1 and 5). In the optically thin limit 
CS column densities are N(CS) ~5xl0 13 cirT 2 and they have 
to be considered as lower limits. Low excitation temperatures 




10 20 30 
upper level energy [K] 



10 20 30 
upper level energy [K] 




10 20 30 40 
upper level energy [K] 
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upper level energy [K] 



Fig. 5. CS rotational-diagrams corrected for line opacity effects 
at each observed position of Fig0] Rotational-diagrams for dif- 
ferent considered CS 7=2-1 line opacities (T2-1) are shown in 
each box. Rotation temperatures for each opacity correction are 
also indicated. 



(T rot ~6-9 K) are also inferred from the rotational-diagrams. 
These values, much lower than expected gas temperatures in a 
PDR, are suggestive of radiative excitation effects in CS lines 
and level populations far from thermalization. Therefore, we 
only use the rotational-diagrams as input for the first iteration 
of a more complex analysis. 

In order to obtain a more detailed overview of the CS ex- 
citation, we have made several statistical equilibrium calcula- 
tions (see Appendix) around the expected physical conditions 
in the Horsehead. In particular, we have run a grid of single- 
component models for Tt=10, 20, 30, 50 and 70 K, n(H 2 )=l0 4 , 
5xl0 4 , 10 s and 5xl0 5 cirT 3 , and^(C5) from 10~ 10 to 10~ 7 . As 
a reference value, the cloud total extinction is assumed to be 
constant and equal to Ay-20 mag in all models, i.e. the spa- 
tial length is changed accordingly. Fig.|6]specifically shows se- 
lected results for Tt=30 K, which gives appropriate absolute 
intensities for the CS lines. In particular, integrated line inten- 
sity ratios of observed lines as a function of CS abundance 
for different densities are shown. Averaged ratios from CS 
single-dish observations are l^f ~0.7, j^~0.2 and |^~0.3. 
Therefore, densities >5xl0 4 cirT 3 are needed to populate the 
CS intermediate-/ levels. On the other hand, for high densi- 
ties (>5xl0 5 cirT 3 ), collisions start to efficiently populate these 
levels and the predicted line ratios involving the CS J -5-4 line 
become much larger than observed. Thus, mean densities are in 
the range n(H 2 ) ^(0.5-1. 0)xl0 5 cirT 3 , i.e. lower than CS crit- 
ical densities (table [5}- Excitation temperatures are predicted 
to be subthermal, T ex <Tt, especially for the highest frequency 
lines. Due to line-trapping, the maximum T ex is reached at the 
center of the cloud, while it graduately drops at both surfaces 
where line photons are optically thin and line trapping is not 
efficient. The only exception is the CS 7=1-0 transition which 
shows an increase of the excitation temperature, T'7 , at both 
surfaces. This rising in T]^ is due to the significant collisional 
excitation coupling from the 7=0 to 7=2 level, and to the large 
radiative de-excitation rates from 7=2 to 7=1 level. At both 
surfaces, where optically thin CS 7=2-1 line photons can eas- 
ily escape from the cloud (T^ 1 decreases), the above excitation 
conditions favor the population of the 7=1 level with respect 
to the 7=0 level. Therefore, T' ~° can reach large suprathermal 
values. A typical example is shown in Fig. IA.2I Thus, within 
this range of parameters and even if physical conditions are ho- 
mogeneous, excitation gradients must be taken into account. 

For these temperatures and densities, T* ^20-30 K and 
n(Hi) ^(0.5-1. 0)xl0 5 cm" 3 , the CS §Ef and f^f line ra- 
tios are better fitted in the interval x(CS) ^(0.2-1.0)xl0~ 8 . 
Nevertheless, the |fy line ratio is systematically predicted to 
be larger than observed in these single-component models. 
Therefore, a more complex density structure and/or additional 
opacity effects in low-/ CS lines may be affecting the observed 
profiles. The latter hypothesis is clearly favored by the fact that 
the C 34 S §5f line intensity ratio is larger (-0.9) than the CS 
|5y ratio (~0.7) and thus closer to the single-component model 
predictions. Since the C 34 S emission is optically thin, radiative 
transfer effects are less important and C 34 S column densities 
can be accurately determined. 
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Table 7. One-component radiative transfer model parameters. Table 8. Two-component radiative transfer model parameters. 



Parameter 



Value 



Parameter 



Value 



n(H 2 ) 

^turb 

X(HCS + ) 



20-25 K 
(7-12)xl0 4 cm" 3 
0.3-0.5 km s" 1 
(3±l)xl0-'° 
(4+2)xl0- n 



C S single-component radiative transfer models for se- 
lected positions within the region have been run (Fig.0. Since 
C 34 S emission can arise from different gas components of 
higher density, not resolved by the large single-dish beam, 
we have modeled each position in spherical geometry. This 
allows to explore different components of different beam fill- 
ing factors. The maximum extinction in the models varies 
from Ay=20 to 12 mag depending on the cloud position. 
These values are consistent with those obtained from single- 
dish 1.2 mm dust continuum emission observations (Teyssier 
et al., 2004, Pety et al. 2005a). Following our previous gen- 
eral excitation calculations we have considered gas temper- 
atures in the range 20-25 K. For these conditions, densities 
between «(//2)=7xl0 4 and 1.2x10 s cirT 3 satisfactorily repro- 
duce the observed C 34 S absolute intensities. Best fits are ob- 
tained for turbulence velocities (see Appendix for the defini- 
tion of \turb) between 0.3 and 0.4 km s (TableQ. Although 
C 34 S is slightly enhanced where C ls O decreases, we have av- 
eraged the 4 positions to find the mean C 34 S abundance in 
the region covered with single-dish observations and found 
;^(C 34 S)=(3+l)xlO~ 10 . Since nucleosynthesis models favor a 
constant galactic 32 S/ 34 S ratio and many observations repro- 
duce the solar ratio within their error bars (Wannier et al. 1980 
; Frerking et al. 1980), especially in local diffuse clouds (Lucas 
& Liszt, 1998), we adopt 32 S/ 34 S=23 here as the isotopic ratio 
in the Horsehead. Therefore, the derived x(C 34 S) abundance 
translates to;^(CS)=(7+3)xlO~ 9 . The same physical conditions 
at each position have been used to model the HCS + 7=2-1 lines 
(see Fig. [8}- Lines are reproduced for an averaged abundance 
of x(HCS + )=(4+2)xlO u , therefore, a CS/HCS+ ^175 abun- 
dance ratio is derived. 

Using the CS abundance inferred from the C 34 S analysis, 
we have now tried to fit the CS lines at each position. Since a 
single-component model does not reproduce the observed line 
ratios and absolute intensities, we have explored other possibil- 
ities. In principle, CS low-/ lines are optically thick and may 
not trace the high density gas revealed by C 34 S, especially if 
the medium is inhomogeneous and dense clumps and a more 
diffuse interclump medium coexist. The same argumentation 
has been used to interpret HCN and H 13 CN observations in the 
Orion Bar PDR (Lis & Schilke, 2003). In addition, it is well 
known that low-/ CS lines may not be a good column den- 
sity tracer if their emission is scattered by a low density halo 
(Gonzalez-Alfonso & Cernicharo 1993). This process can be a 
common effect in optically thick lines of high-dipole moment 
molecules such as CS or HCO + (Cernicharo & Guelin 1987). 
In this scenario, the CS J-3-2 and 2-1 lines from the dense 



I* 

n(H 2 ) 
dense component 
(filling factor) 

dense component 
X(CS) 
S/H 



20-25 K 
(3-7)xl0 4 cnr 3 
(2-6)xl0 5 cm" 3 
0.3 

0.3-0.4 km s- [ 
0.2-0.3 km s- 1 

(7±3)xl0- 9 
(3.5±1.5)xlQ- 6 



medium will be attenuated and scattered over larger areas than 
the true spatial extend of the dense clumps. This possibility has 
been analyzed in more detail in the next section. Fortunately, 
observations of the CS /=5-4 line allow to directly trace the 
dense clumps more safely (Table In particular, we found 
that these lines can only be reproduced with denser gas com- 
ponents, n(H 2 )= (4±2)xl0 5 cm 4 , not resolved by the ~10" 
beam of the IRAM-30m telescope at -250 GHz. Note that the 
CS /=5^1 line widths are fitted if the turbulent velocity in the 
denser gas is ~0.2 km s _1 , a factor 2 lower than the one required 
by the C 34 S /=3-2 and 2-1 lines (Fig.0. Thus, a different spa- 
tial origin for this line emission is favored. 

At this stage we have a general knowledge of the CS and 
HCS + excitation and abundance in the region. In the following 
sections we concentrate in the photochemistry of these species. 
Only higher angular observations provide the appropriate linear 
scale to resolve the most important physical gradients in the 
PDR edge. Hence, interferometric observations should allow a 
better comparison with chemical predictions. 

5.2. The PDR edge 

PdBI C 12 7=2-1, 1-0, C I8 7=2-1, and CS /=2-l observa- 
tions along the direction of the exciting star (at 6y- 0") are 
shown in Fig. [3] Here we take these spectra as representa- 
tive of the PDR edge and try to constrain its physical condi- 
tions through a combined analysis of photochemical and radia- 
tive transfer models. Both models use a unidimensional plane- 
parallel description of the geometry. Although some physical 
processes require more complex geometries, the main physi- 
cal and chemical gradients across the illuminated direction can 
be consistently described in this way. Plane-parallel geometry 
was judged to be the best approach for this edge-on PDR since 
H2 and PAH emissions are only observed at the illuminated 
edge and not deeper inside the cloud (Habart et al. 2005). 

In this analysis, we have used the PdBI CS 7=2-1 and 
C ls O 7=2-1 lines. As low-/ 12 CO optical depths are very high, 
they do not trace the bulk of material. The intensity peak of 
these lines only provide a good estimation of the CO excitation 
temperatures (i.e. a lower limit to the gas temperature). Since 
the asymptotic brightness temperature of CO /=l-0 lines is 
~30 K, we take this value as the minimum of Tj. in the PDR. 
We note that lower temperatures do not reproduce the observed 
line intensities. For the rest of the (warmer) positions closer to 
the PDR edge, the gas temperature was determined by solv- 
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ing the thermal balance. The predicted gas temperature in the 
density peak is ~50 K while it rises to ~200-250 K in the 
H2 emitting regions where the density is n H ^10 3 -10 4 cm" 3 . 
More exact temperature values require observations of higher- 
7 CO lines at comparable spatial resolution. We are currently 
analysing 13 CO 7=3-2 data from the SMA interferometer. 

Regarding the density structure, both the observed H2 and 
PAH mid-IR emission, together with their spatial segregation, 
are much better reproduced with a steep density gradient than 
with an uniform density (Habart et al. 2005). The same den- 
sity gradient is needed to correctly reproduce the observed off- 
set between the small hydrocarbons (Pety et al. 2005a) and H2 
emission (where the density is not at its peak). Therefore, in 
order to reproduce PdBI observations of CS and C 18 0, a steep 
power-law density gradient at the illuminated regions and a 
step-density in the more shielded region have been assumed. 

The following methodology was carried out: a full PDR 
model with Horsehead standard conditions (see section 14. It 
was run with a particular choice of the density gradient de- 
scribed in Eq. |2] Afterwards, the PDR output was used as in- 
put for the nonlocal radiative transfer calculation in a fash- 
ion described in appendix IA.2I In this way, physical param- 
eters can be tuned more accurately by iteration of differ- 
ent radiative transfer models. Once better parameters have 
been found, a new PDR computation is performed with this 
choice of physical parameters. Hence, the most appropri- 
ate physical and chemical description of the PDR edge was 
found through the PDR model^transfer model^check with 
observations^transfer model^PDR model iterative process. 
Therefore, synthetic CS and C 18 abundance profiles are con- 
sistently computed as a function of the edge distance Sx (in 
arcsec) and directly compared with observations. 

Different PDR spatial depths were investigated. Depending 
on the adopted density profile, the spatial depth Idepth is deter- 
mined by the line of sight visual extinction. However, the Ay 
value depends on the method used to measure it. If optically 
thin 1.2 mm dust emission is used (Teyssier et al. 2004; Pety 
et al. 2005a; Habart et al. 2005), the resulting column densi- 
ties depend on the usually unknown grain properties and on 
the assumed temperature. Taking into account our poor knowl- 
edge of the cloud thermal structure, a factor ~2 of uncertainty 
in Ay can be assumed. In addition, the angular resolution of 
millimeter continuum observations is at least a factor ~2 worse 
than PdBI molecular line observations. Due to the steep de- 
crease of the density towards the edge, and due to the ~11" 
beam of 1.2 mm continuum observations, the observed emis- 
sion peak will appear deeper inside the cloud, shifted a few 
arcsec from the real density peak (which is closer to the edge). 
Therefore, together with the PDR edge location, the exact peak 
density position can also be uncertain by a few arcsec. Finally, 
beam dilution has to be also taken into account when compar- 
ing single-dish versus interferometric observations. Here we 
have chosen Idepth =0.05-0.1 pc, which implies extinction peaks 
around Ay —15-30 magnitudes. These values are expected in 
compact globules (Reipurth & Bouchet 1984). Since CS and 
C ls O excitation and line transfer are quite different, the fol- 
lowing combined analysis provides an accurate description of 



the edge density structure. The empirical density profile in the 
models, nH-n(H)+2n(Hi), as a function of Sx is: 

' «h(0) + [n H (8x{) - n H (0)] (^f ; 6x1 > Sx > 
n H (Sx) = - n H (8x\); 6x 2 > Sx > Sxi ( 2 ) 

«h(c«2); Sx > 5x2 
where 6x is the distance away from the PDR edge, «h(0) is 
the ambient density at the edge, and nniSxi) and naiSxz) are 
constant densities in the 6x2 > Sx > Sx\ and Sx > 6x2 re- 
gions respectively. Selected photochemical models are shown 
in Fig. [9] The normalized population of the H2 v=l, 7=3 level 
is shown in the upper panel and is used to place the £x-axis 
origin of the models and thus to accurately check with obser- 
vations. Although some uncertainty in the location the PDR 
edge exists, we place the peak of this curve at the maximum of 
observed H2 1-0 S(l) 2.12 /im excited line (Sx ~10"; Habart 
et al. 2005). Best models are found for a peak density around 
«h(^i)-2x10 5 crrT 3 . This density is reached in a length of 
~2.5"-5" (or 5-10xl0~ 3 pc) and stays constant in a length of 
SX2 - Sx\ <20" (or 0.04 pc). In order to fit the smooth de- 
crease of C ls O emission and also of the 1.2 mm continuum 
emission, the density has to decrease again by at least a factor 
~2. We have simply modeled this as a step-function for Sx>Sx2 
and decrease the density to «h((5x2)=10 5 crrT 3 . We have cho- 
sen Sx\ = 12" and 8x2 = 30". Our models confirm that high 
density gas and a large gradient slope, [3 ~3-4, are needed to 
reproduce the PdBI and H2 observations (Habart et al. 2005), 
although we found a slightly smaller gradient scale length. 

As proposed by Habart et al. (2005) the PDR edge can be 
slightly inclined with respect to the line of sight by a small an- 
gle ip. In plane-parallel geometry, the maximum inclination can 
be estimated assuming that the observed spatial extend of the 
H2 emission, dH 2 , is mainly due to the projection of Idepth in the 
plane of the sky, thus sinip ^ dH 2 l Idepth- Since du 2 ~0.01 pc, 
an inclination angle <p ~5°, has been considered in the radia- 
tive transfer models (see Appendix IA.2> . As expected, even 
such a small inclination shifts the emission peak significantly 
and should therefore be taken into account. Fig. ^3 shows the 
PdBI C 18 line observations and the combined PDR+transfer 
modeling including such geometrical effects. The agreement is 
excellent, probably favored by the well-established CO photo- 
chemistry (Fig. [9} and because C I8 7=2-1 lines do not show 
complex radiative transfer effects (j2-\ ~0.8). 

To analyse the spatial distribution of the CS abundance pre- 
dicted by photochemical models at the PDR edge we have also 
tried to fit the PdBI CS 7=2-1 lines. Fig.|5]shows the effects 
of different sulfur abundances; S/H=2xl0~ 5 and S/H=2xl0~ 6 . 
Fig. \^\ (no inclination) and Fig. ^1 (inclination considered) 
show the resulting synthetic CS map, using S/H=2xl0~ 6 and 
a minimum gas temperature of 30 K, over PdBI observations 
at two constant Sy cuts (Sy=30" and 0"). Contrary to C ls O, the 
CS emission detected with the PdBI at a fixed Sx near the edge 
shows an emission gradient in the Sy direction, e.g. line peaks 
are brighter as Sy increases. As a consequence model predic- 
tions fit better the Sy-30" cut than the Sy-Q" one. Besides, 
larger gas phase sulfur abundances are obtained if the bulk of 
the gas in the PDR edge is warmer, i.e. minimum gas tempera- 
tures of ~100 K (Fig. [9] right panel). This may be an indication 
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Fig. 13. CS and C S synthetic line profiles for a cloud 
with a depth of 0.1 pc, T t =30 K, n(H 2 )=10 5 cirT 3 and 
;tf(CS)=7xlO~ 9 (thick curves). Thin curves show the result- 
ing spectra if the same cloud is surrounded by different dif- 
fuse halos (3: «(H 2 )=5xl0 3 cm- 3 , 2: n(H 2 )=8xl0 3 cirT 3 and 1: 
«(H 2 )=lxl0 4 cm 3 ). The CS abundance in the cloud is deter- 
mined more precisely from CS high-/ and C 34 S low- J obser- 
vations, otherwise it is underestimated. Note that the intensity 
levels are comparable to those observed in the Horsehead. 

of larger temperatures at the cloud edge and lower sulfur deple- 
tions. Note that an accurate estimation of the CS abundance at 
high resolution from a single PdBI line is not straightforward. 
Such determination requires aperture synthesis observations of 
additional CS lines to have a minimum idea of the CS excita- 
tion in different positions. 

In addition, since C 34 S observations at the same high- 
angular resolution were not available, we could not estimate ad- 
ditional opacity effects in previous PdBI CS models (t 2 _i >2). 
In the following we have tried to estimate the worse possible 
scenario affecting the CS lines in the line-of-sight, i.e. the 
presence of a surrounding low density halo. Of course, the 
greatest effect can appear in the shielded regions where the 
gas column density is largest. Therefore we modelled a typi- 
cal position where CS is well spatially resolved with the fol- 
lowing parameters: ld ep rh-0-l pc, Tt=30 K, n(H 2 )=10 5 cm -3 , 
v„, ri =0.35 km s" 1 and^(CS)=7xl(T 9 , (the averaged CS abun- 
dance obtained from the detailed CS and C 34 S excitation anal- 
ysis of previous section). We consider in addition that a low 
density halo of diffuse gas with the same ;f (CS) surrounds the 
region. We take Ta = 10 K, v (l „-/,=0.7 km s _1 and densities in the 
interval n(H 2 )=(5-10)xl0 3 cirT 3 . The same modeling was car- 
ried out for C 34 S. Fig. 1131 shows model results. As expected, 
a low density halo efficiently self-absorbs CS line photons in 
the most opaque lines, i.e. the low-/ CS lines. As a result, the 
observed CS line intensities are attenuated and abundances can 
be easily underestimated. However, this effect can be differ- 
ent at different positions, since the line opacity also changes. 
Apart from uncertainties in sulfur chemistry or instrumental ef- 



fects in interferometric observations, diffuse gas can also con- 
tribute to explain differences between models and observations 
in Figs.^2 an dEl Since, optically thick lines are affected by 
this effect (Figll3l>. only the observation of 13 CS or C 34 S iso- 
topologues can help to provide more accurate abundance de- 
terminations. In the following, the CS chemistry is analyzed in 
more detail. 

5.3. CS chemistry and S-abundance 

Predicted C ls O/CO/C/C + and CS/HCS + /S/S + structures for a 
unidimensional PDR with Horsehead standard conditions are 
shown in Fig. |5] (see section |4~TV Variation of the sulfur ele- 
mental abundance almost does not affect the CO or C I8 abun- 
dance profiles, but it slightly modifies the predicted C/C + abun- 
dance profiles because charge transfer reactions between C + 
and S, and between C and S + are clearly altered by the sul- 
fur depletion. The following results are of course determined 
by our present knowledge and uncertainties on S-chemistry, 
and on reaction rates at different temperatures. According to 
the latest ion storage ring experiments (Montaigne et al. 2005), 
only 19% of the HCS + dissociative recombination results in 
CS + H while the fracture of the C-S bond dominates the 
dissociation (81%). Since these experiments can not sepa- 
rate the contribution of the CH + S or SH + C channels in 
the latter process, we have adopted the same branching ra- 
tio (0.405) for both channels. The reaction rate coefficient is 
£ M (i/CS + ) /flJf =9.7xl0- 7 (T/300r - 57 cm 3 s _1 . We have also 
included the latest OCS + dissociative recombination rates from 
Montaigne et al. (2005). The CS + O production channel now 
occurs at a rate 3 times slower than in previous chemical net- 
works. All these modifications clearly influence the amount 
of CS formed from a given sulfur abundance, and thus the 
sulfur depletion estimations. Fig. [5] (left panel) also shows 
the effect of adopting the older HCS + and OCS + dissocia- 
tive recombination rates and branching ratios. In particular, 
fczj fi (#CS + ) jW =5.8xl0- 8 (T/300r a75 cm 3 s -1 . However, since 
HCS + + <?~ — > CS + H was the only channel considered and 
the OCS + + e~ — > CS + O process was faster, smaller sulfur 
abundances were required to obtain the same CS abundances. 

In the most external layers of the cloud, still dominated by 
the FUV radiation field, CS is predominantly formed by HCS + 
dissociative recombination and principally destroyed by pho- 
todissociation and charge transfer with H + . Once the gas is 
shielded, OCS + dissociative recombination and reaction of C 
with SO also contributes to CS formation, while its destruc- 
tion is now governed by ion-molecule reactions, mainly with 
HCO + but also with H30 + . These last two reactions with abun- 
dant molecular ions return HCS + again. The peak abundance 
of HCS + occur at Ay <2 mag, where it is formed by reaction 
of CS + with H 2 and destroyed by dissociative recombination. 
For this reason, an order of magnitude change in Icdr(HCS + ) 
clearly modifies its peak abundance in the outer PDR layers. 
In the more shielded regions, HCS + destruction is dominated 
by dissociative recombination and reaction with atomic oxygen 
to form HCO + and OCS + . Since the predicted CS abundance 
scales with S/H, and CS formation is dominated by HCS + dis- 
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sociative recombination, we have used our CS/C 34 S/HCS + ob- 
servations and modeling to estimate S/H. 

Fig. ^] shows results of a grid of photochemical models 
for different sulfur elemental abundances from S/H=1(T 8 to 
2xl0~ 5 , using the latest HCS + and OCS + dissociative recom- 
bination rates. CS and HCS + abundances with respect to H2 
are shown as a function of S/H at two different PDR posi- 
tions (A v ~10 and ~2 mag respectively; see Fig.|5}. Densities 
at these positions are the same, n{H2)=10 5 cm -3 , but we have 
taken different PDR positions in order to plot the HCS + maxi- 
mum abundance and to get the CS/HCS + ratio closer to obser- 
vations. Inside the cloud, the predicted maximum HCS + abun- 
dances are a factor ~3 lower than observed. Horizontal shaded 
regions mark the CS and HCS + abundances derived from ob- 
servations and radiative transfer modeling. For clarity, HCS + 
abundances have been multiplied by a factor of 1000. Finally, 
the vertical shaded region shows the estimated sulfur elemental 
abundance in the Horsehead derived from the overlap region 
between observed and predicted abundances. We derive S/H 
~(3.5±1.5)xl0~ 6 as the mean value for the PDR. Note that 
CS is used for the upper limit and HCS + for the lower limit. 
However, according to the inferred HCS + abundance, larger 
sulfur abundances are still possible. 




10 8 10 7 10 6 10 5 

gas phase S/H 

Fig. 14. Photochemical model predictions for the physical and 
FUV illuminating conditions prevailing in the Horsehead PDR 
showing the CS and HCS + abundance as a function of the sul- 
fur gas phase abundance. Horizontal shaded regions show the 
CS and HCS + abundances derived from the single-dish obser- 
vations and radiative transfer modeling. Note that for clarity 
HCS + abundances have been multiplied by a factor of 1000. 
The shaded vertical region shows the estimated sulfur abun- 
dance in the Horsehead nebula derived from the constrained 
fits of CS and HCS + abundances. 



6. Discussion 

Our multi-transition single-dish and aperture synthesis obser- 
vations and modeling of CS and related species allow us to 
constrain the sulfur gas phase chemistry in the Horsehead PDR 
and it also gives some insights on the dense gas properties. 



6.1. Densities 

The densities found in this work, n{Hi) ^10 5 cirT 3 , are larger 
to those inferred from previous studies based on single-dish 
CO observations (Abergel et al. 2003; Teyssier et al. 2004). 
This may be the indication of an inhomogeneous medium char- 
acterized by a interclump medium (well traced by CO) and a 
denser clump medium (better traced by high dipole molecules). 
Both high densities and inhomogeneous medium are common 
in other PDRs such as the Orion Bar (Lis & Schilke 2003). 
In particular, we have shown that unresolved gas components 
up to n(Hj) ^(2-6)xl0 5 cirT 3 are required to explain the CS 
J-'S-A line emission in the Horsehead. However, Abergel et 
al. (2003) did not find inhomogeneities in analysing ISOCAM 
images of the Horsehead. Nevertheless, they noted that dumpi- 
ness at scales smaller than the upper limit of the FUV penetra- 
tion depth (~0.01 pc) could not be excluded. Our best mod- 
els of the CS /=5— 4 line emission require an unresolved com- 
ponent with a radius of ~5xl0~ 3 pc. This component can of 
course be further fragmented itself. Nevertheless, it is difficult 
to distinguish between dumpiness at scales below ~0.01 pc and 
the presence of a lower density envelope surrounding the cloud. 
Since CO 7=1-0 and 2-1 line opacities easily reach large val- 
ues, their observed profiles are formed in the very outer lay- 
ers of the cloud and thus they can arise from the most dif- 
fuse gas (n{Hi) ~5xl0 3 cm -3 ). Interferometric observations 
of intermediate-/ lines of high dipole species such as CS or 
HCO + will help to clarify the scenario. 

The high angular resolution provided by PdBI CS and C I8 
observations reveals that the Horsehead PDR edge is character- 
ized by steep density gradient rising from ambient densities to 
n(H-i) ~10 5 cirT 3 in a length of ~0.01 pc and kept roughly con- 
stant up to ~0.05 pc, where the density decreases again at least 
a by factor 2. The exact density values still depend on the as- 
sumed cloud depth and temperatures. In any case, the inferred 
shell of dense molecular gas has high thermal pressures ~(5- 
10)xl0 6 K cirT 3 and this can be the signature of the processes 
driving the slow expansion of the PDR. Therefore, the most 
shielded clumps undergo effective line cooling and the regions 
of lower density should be compressed due to their lower inter- 
nal pressure. Recent hydrodynamical simulations of the expan- 
sion of ionization and dissociation fronts around massive stars 
also predict that a high density molecular shell (10-100 times 
the ambient density) will be swept up behind the ionization 
front (Hosokawa & Inutsuka 2005 a&b). The density, pressure 
and temperature profiles and values predicted by these simula- 
tions at ~0.5 Myr (the Horsehead formation timescale derived 
from its velocity gradients by Pound et al. 2003 and Hily-Blant 
et al. 2005) qualitatively reproduce the values inferred from 
our molecular line observations and modeling. Hence, a shock 
front driven by the expansion of the ionized gas is probably 
compressing the cloud edge to the high densities observation- 
ally inferred in this work. Specific hydrodynamical simulations 
for the particular source physical conditions and comparison 
with observations will be appreciated. As noted by Hosokawa 
& Inutsuka (2005), the dynamical expansion of a HII region, 
PDR and molecular shell in a cloud with a density gradient 
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has not been studied well. We suggest the Horsehead PDR as a 
good target. 



6.2. Temperatures 



Molecular excitation, radiative transfer and chemical mod- 
els are used to derive realistic abundances. The gas tempera- 
ture impacts many aspects of these computations (e.g. chem- 
ical reaction rates and collisional excitation), and thus, the 
density and abundance uncertainties also reflects our incom- 
plete knowledge of the thermal structure. The problem is not 
straightforward, since a steep temperature gradient is also ex- 
pected in PDRs, and also because the most appropriate tracers 
of the warm gas lie at higher frequencies. The Horsehead PDR 
may not be an extreme case, since its FUV radiation field is 
not very high and photoelectric heating alone will not heat the 
gas to high temperatures as long as the gas is FUV-shielded. 
Nevertheless, our thermal balance calculations quickly lead to 
Tt ^10 K. According to observations, this temperature is too 
low, especially in the first 6x ~30" representing the PDR edge. 
In this work we have (observationally) adopted Tt=30 K as the 
minimum gas temperature in our PDR calculations. In forth- 
coming works we will concentrate on the thermal structure of 
the PDR. Here we only note that either the cooling is not so ef- 
fective and/or extra heating mechanisms need to be considered. 
The cosmic ray ionization rate was also increased by a factor 
~5 but it only modifies the gas temperature by AT^ ~4 K. 

Under these circumstances we have to conclude that the 
gas, or at least a fraction of it, is likely to be warmer than pre- 
dicted. We note that this problem is not new. Again, high-/ 
13 CO and NH3 observations have previously probed the exis- 
tence of warm gas (~150 K) in regions where standard heating 
mechanisms fail to predict those values (Graf et al. 1990; Batrla 
& Wilson 2003). More recently, Falgarone et al. (2005) re- 
ported ISO observations of H2 in the lowest five pure rotational 
lines S(4) to S(0) (8 fim to 28 /an) toward diffuse ISM gas. The 
observed S(1)/S(0) and S(2)/S(0) line ratios are too large to 
be compatible with the PDR models. These authors suggested 
that MHD shocks (Flower & Pineau des Forets 1998) or mag- 
netized vortices, which are natural dissipative structures of the 
intermittent dissipation of turbulence (Joulain et al 1998), lo- 
cally heat the gas at temperatures up to * 1000 K. These struc- 
tures add two heating sources: i) viscous heating through large 
velocity shear localized in tiny regions and ii) ion-neutral drift 
heating due to the presence of magnetic fields (ambipolar dif- 
fusion). As shown by Falgarone et al (2006), these dissipative 
structures are also able to trigger a hot chemistry, that is not 
accessible to models that do not take into account the gas dy- 
namics. These results suggest that additional mechanical heat- 
ing processes are at work. We propose that the shock waves 
driven by the expansion of the HII region and PDR compress 
the molecular gas in the cloud edge and provide it with an ad- 
ditional heating source. 



6.3. CS and HCS + chemistry 

According to the last (but not least) molecular data affect- 
ing CS chemistry and excitation, the mean CS abundance in 
the Horsehead PDR, ;KCS)=(7±3)xlO~ 9 , implies a gas sulfur 
abundance of S/H ~(3.5+1.5)xl0~ 6 , only a factor <4 smaller 
than the solar value ( Asplund et al. 2005). Even lower sulfur 
depletion values are possible if the gas is significantly warmer 
than 30 K. Thus, the gas phase sulfur abundance is very close to 
the undepleted value observed in the diffuse ISM and not to the 
depleted value invoked in dense molecular clouds (e.g. Millar 
& Herbst 1990). However, the observed CS/HCS+ ^175 abun- 
dance ratio can only be reproduced by photochemical models 
by considering the HCS + peak abundance, otherwise, larger ra- 
tios (~10 3 ) are predicted. Therefore, either the observed HCS + 
only traces the surface of the cloud where its abundance peaks, 
or chemical models underestimate the HCS + production rate. 
In any case, the predicted CS/HCS + abundance ratio scales 
with the gas phase sulfur abundance. The largest ratios are ex- 
pected at the lowest sulfur depletions. However, the observed 
CS/HCS+ -10 ratio in the diffuse ISM (Lucas & Listz 2002) is 
even lower than in the Horsehead. Thus, we have to conclude 
that present chemical models still fail to reproduce the observed 
CS/HCS + abundance ratio, at least in the stationary regime. 
Time-dependent photochemical computations may also help to 
understand the dynamical expansion of the dissociation front 
and the evolving molecular abundances. Besides, Gerin et al. 
(1997) noted that larger HCS + abundances are expected if the 
gas is in a high ionization phase. We have computed that if the 
cosmic ray ionization rate is increased by a factor of ~5, the 
predicted HCS + abundance inside the cloud (Ay=10 mag) in- 
terestingly matches our inferred value and the CS/HCS + abun- 
dance ratio gets much closer to the observed ratio without the 
need of taking the HCS + abundance peak. 

For the physical and FUV illuminating conditions prevail- 
ing in the Horsehead PDR, most of the gas phase sulfur is 
locked in S + for Ay <2 mag and x(HCS + ) reaches its maxi- 
mum value. Besides, the derived gas phase sulfur abundance 
is large enough to keep ;^(e~)>10~ 7 for Ay<3.5 mag. HCS + 
and S 11 recombination lines trace the skin of molecular clouds 
where S + is still the dominant form of sulfur. In the scenario 
proposed by Ruffle et al. (1999), these S + layers will be re- 
sponsible of the sulfur depletion due to more efficient sticking 
collisions on negatively charged dust grains than in the case 
of neutral atoms such as oxygen. Even in these regions, still 
dominated by photodissociation, CS and HCS + abundances are 
quickly enhanced compared to other sulfur molecules. In fact, 
we predict that CS is the most abundant S-bearing molecule in 
the external layers where S + is still more abundant than neutral 
sulfur. These results are consistent with our PdBI detection of 
CS close to the PDR edge and show that CS is a PDR tracer. 
These findings are consistent with observations of S-bearing 
species in the diffuse ISM where CS is more abundant than 
S0 2 , H 2 S and SO (Lucas & Listz 2002). 

Between Ay ~2 and ~8 mag the S + abundance smoothly 
decreases. Since S + is a good source of electrons, the electronic 
fractionation also decreases accordingly. HCS + , and thus CS, 
present an abundance minimum in these layers. Neutral atomic 
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sulfur is now the most abundant S-bearing species. Therefore, 
observations of the [S i]25 //m fine structure line will basi- 
cally trace these intermediate layers of gas where S-bearing 
molecules have not reached their abundance peak. However, 
the exctitation energy of the [S i]25 //m line (the upper level 
energy is ~570 K) is too high compared to the thermal en- 
ergy available in the regions where the neutral sulfur abun- 
dance peaks (7* ^30 K) and no detectable emission is expected. 
In fact, no Spitzer/IRS line detection has been reported in the 
Horsehead (L. Verstraete, private com.). However, since most 
of the neutral atomic sulfur will remain in the ground-state, 
the presence of a background IR source (e.g. in face-on PDRs) 
may allow, with enough spectral resolution and continuum sen- 
sitivity, the detection of the [S i]25 fim line in absorption. 

On the other hand, sulfur in diffuse ionized gas outside the 
molecular cloud is in the form of sulfur ions. Mid-IR [S m] 
fine structure lines have been detected around the Horsehead 
with IRS/Spitzer (L. Verstraete, private com.). In the shielded 
gas, sulfur is mostly locked in S-bearing molecules together 
with a smaller fraction in atomic form. Our models predict 
that species such as SO will be particularly abundant. Jansen 
et al. (1995) also noted that the low gas phase sulfur abundance 
needed to explain the CS abundance in the Orion Bar PDR was 
incompatible with the observed high H2S/CS~0.5 abundance 
ratio. Therefore, a complete understanding of the sulfur chem- 
istry will only be achieved when all the major sulfur molecules 
can be explained. In a forthcoming paper we analyse the photo- 
chemistry, excitation and radiative transfer of several S-bearing 
molecules detected by us in the Horsehead PDR. 

7. Summary and Conclusion 

We have presented interferometric maps of the Horsehead PDR 
in theCS 7=2-1 line at a 3.65"x3. 34" resolution together with 
single-dish observations of several rotational lines of CS, C 34 S 
and HCS + . We have studied the CS photochemistry, excitation 
and radiative transfer using the latest HCS + and OCS + disso- 
ciative recombination rates (Montaigne et al. 2005) and CS col- 
lisional cross-sections (Lique et al. 2006). The main conclu- 
sions of this work are as follows: 

1 . CS and C 34 S rotational line emission reveals mean densities 
around n(// 2 )=(0.5-1.0)xl0 5 cm" 3 . The CS /=5^ lines 
show narrower line widths than the low-/ CS lines and 
require higher density gas components, ~(2-6)xl0 5 cirT 3 , 
not resolved by a ~10" beam. These values are larger than 
previous estimates based on CO observations. It is likely 
that dumpiness at scales below ~0.01 pc and/or a low den- 
sity envelope play a role in the CS line profile formation. 

2. Nonlocal, non-LTE radiative transfer models of opti- 
cally thick CS lines and optically thin C 34 S lines pro- 
vide an accurate determination of the CS abundance, 
^(C5)=(7±3)xl0 -9 . We show that radiative transfer and 
opacity effects play a role in the resulting CS line pro- 
files but not in C 34 S lines. Assuming the same phys- 
ical conditions for the HCS + molecular ion, we find 
X(HCS + )=(4±2)xl0- n . 



3. According to photochemical models, the gas phase sul- 
fur abundance required to reproduce these CS and HCS + 
abundances is S/H=(3.5+1.5)xl0~ 6 , only a factor ~4 less 
abundant than the solar elemental abundance. Larger sulfur 
abundances are possible if the gas is significantly warmer. 
Thus, the sulfur abundance in the PDR is very close to 
the undepleted value observed in the diffuse ISM. The pre- 
dicted CS/HCS + abundance ratio is close to the observed 
value of ~ 175, especially if predicted HCS + peak abun- 
dances are considered. If not, the HCS + production is un- 
derestimated unless the gas is in a higher ionization phase, 
e.g. if the cosmic ray ionization rate is increased by ~5. 

4. High angular resolution PdBI maps reveal that the CS emis- 
sion does not follow the same morphology shown by the 
small hydrocarbons emission the PDR edge. In combina- 
tion with previous PdBI C ls O observations we have mod- 
eled the PDR edge and confirmed that a steep density gra- 
dient is needed to reproduce CS and C ls O observations. 
The resulting density profile qualitatively agrees to that pre- 
dicted in numerical simulations of a shock front compress- 
ing the PDR edge to high densities, n(H2)-10 5 cirT 3 , and 
high thermal pressures, ^(5-10)xl0 6 K cirT 3 . 

5. Conventional PDR heating and cooling mechanisms fail to 
reproduce the temperature of the warm gas observed in the 
region by at least a factor ~2. Additional mechanical heat- 
ing mechanisms associated with the gas dynamics may be 
needed to account for the warm gas. The thermal structure 
of the PDR edge is still not fully constrained from obser- 
vations. This fact adds uncertainty to the abundances pre- 
dicted by photochemical models. 

We have shown that many physical and chemical variations 
in the PDR edge occur at small angular scales. In addition, the 
molecular inventory as a function of the distance from the il- 
luminating source can only be obtained from millimeter inter- 
ferometric observations. High angular resolution observations 
contain detailed information about density, temperature, abun- 
dance and structure of the cloud, but only detailed radiative 
transfer and photochemical models for each given source are 
able to extract the information. A minimum description of the 
source geometry is usually needed. Future observations with 
ALMA will allow us to characterize in much more details many 
energetic surfaces such as PDRs. 
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Fig. 6. Grid of CS single-component models assuming T^= 30 K and a fixed extinction of 20 mag. Panels show different line 
ratios as a function of x(CS). Each panel correspond to a single density, from n(H2)=10 4 to 5xl0 5 cirT 3 . Mean observed ratios 
are §Ef-0.7, §5y^0.2 and ^0.3. 
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Fig. 7. Radiative transfer models for CS and C 34 S discussed in the text (curves) that best fits the IRAM-30m observations 
(histograms). Offsets in arcsec refer to the (0,0) position of the C 18 0(2-l) map (see Fig|4}. Predicted line profiles have been 
convolved with the telescope angular resolution at each frequency. Intensity scale is in main beam temperature. 
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Fig. 8. IRAM-30m HCS + (2-l) single-dish observations (histograms) at different positions of the Horsehead. Offsets in arcsec 
refer to the (0,0) position of the C I8 0(2-1) map (see Fig|4} Radiative transfer models for HCS + at selected positions are also 
shown (curves). Predicted line profiles have been convolved with the telescope angular resolution at each frequency. Intensity 
scale is in main beam temperature. 
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Fig. 9. Photochemical models using a unidimensional PDR code for two different sulfur gas phase abundances, (S/H=2xl(T 5 ; 
left) and the minimum value found for the Horsehead (S/H=2xl0~ 6 ; right). The predicted normalized population of the H2 v=l, 
J-3 level is shown in the upper panel and is used to place the £x-axis origin for the models. The peak of this curve is placed 
at the maximum of the observed H2 1-0 S(l) 2.12 pm line emission (5x~10"; Habart et al. 2005). Next panel shows the density 
profile (n H - n(H) + 2n(H 2 ) in cm 3 ) used in the PDR calculations that better fits the CS and C ls O IRAM-PdBI observations. 
Next panel shows the gas temperature (in K) consistently computed in thermal balance until reaches a minimum value of 30 K. 
Lower panels show the spatial variation of C 18 0/CO/C/C + and CS/HCS + /S/S + abundances (relative to n#) across the PDR. The 
far-UV radiation field is ;f=60 times the Draine field. Chemical rates are those of the Ohio State University (osu) gas-phase 
chemical network (September 2005 release) plus several modifications (see text). Bottom left panel shows the effect of using 
the older rate and branching ratios for the HCS + dissociative recombination on the CS and HCS + abundances (dashed curves). 
Bottom right panel shows the effect of using a minimum gas temperature of 100 K in the chemistry. Lower CS abundances and 
thus larger S/H values are possible as the temperature increases (dashed curves). 
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Fig. 10. IRAM-PdBI C ls O 7=2-1 spectra along the direction of the exciting star at 6y = 0" (histograms). Radiative transfer 
models using the output of PDR models for C ls O (blue curve) for a density gradient and physical conditions discussed in the 
text. Lower panel shows inclination effects assuming that the PDR is inclined relative to the line of sight by a <p=5° angle. 
Modeled line profiles have been convolved with an appropriate gaussian beam corresponding to each synthesized beam. Intensity 
scale is in brightness temperature and abscissa in LSR velocity. 
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Fig. 11. IRAM-PdBI CS 7=2-1 spectra along the direction of the exciting star at 5y - 30" (upperpanel) and 6y = 0" 
(lowerpanel). Radiative transfer models using the output of PDR models for CS (red curve) for a density gradient and phys- 
ical conditions discussed in the text (assuming that the PDR is not inclined relative to the line of sight). Modeled line profiles 
have been convolved with an appropriate gaussian beam corresponding to each synthesized beam. Intensity scale is in brightness 
temperature and abscissa in LSR velocity. 
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Fig. 12. IRAM-PdBI CS 7=2-1 spectra along the direction of the exciting star at 5y - 30" (upperpanel) and 6y = 0" 
(lowerpanel). Radiative transfer models using the output of PDR models for CS (red curve) for a density gradient and phys- 
ical conditions discussed in the text (assuming that the PDR is inclined relative to the line of sight by a ip=5° angle). Modeled 
line profiles have been convolved with an appropriate gaussian beam corresponding to each synthesized beam. Intensity scale is 
in brightness temperature and abscissa in LSR velocity. 
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Fig. A.l. Plane-parallel geometry for a cloud isotropically illu- 
minated by the cosmic microwave background at both surfaces. 

Appendix A: Nonlocal, non-LTE radiative transfer 

Radiative transfer in a medium dominated by gas phase 
molecules and dust grains requires the solution of the radiative 
transport (RT) equation for the radiation field together with the 
equations governing the relative level populations of the con- 
sidered species. In the case of rotational line emission (far-IR 
to mm domain), scattering from dust grains can be usually ne- 
glected from the RT equation and steady state statistical equi- 
librium can be assumed for molecular populations. However, 
physical conditions in ISM clouds are such that molecular ex- 
citation is usually far from LTE. Therefore, a minimum treat- 
ment of the nonlocal coupling between line+continuum radia- 
tion and level populations is required. In this appendix we de- 
scribe in more detail the simple model developed for this work. 

A.1. Monte Carlo methodology for gas and dust 
radiative transfer in plane-parallel geometry 

The Monte Carlo methodology or its modifications is a simple 
and widely adopted approach when one has to deal with mod- 
erately thick lines and flexibility to explore different geome- 
tries is required (see van Zadelhoff et al. 2002 and references 
therein). In this work, the classical description of the Monte 
Carlo approach for non-LTE line transfer (Bernes 1979) has 
been extended to include the dust emission/absorption and their 
effect on the source function. The code was originally devel- 
oped in fortranW for spherical symmetry (Goicoechea 2003) 
and has been enlarged to semi-infinite plane-parallel geome- 
try (from face- to edge-on). Thus, numerical discretization is 
transformed from spherical shells to slabs. The model includes 
illumination from the cosmic background at both surfaces. 

The variation of the radiation field intensity along any pho- 
ton path s is related to the emission and absorption properties 
of the medium (scattering neglected) through 

-j L =j v -a v I v (A.l) 
as 

where a v [cm -1 ] and j v [erg s _1 crrT 3 Hz T 1 sr _1 ] are the to- 
tal (gas+dust) absorption and emissivity coefficients at a given 
frequency v. The normal path to any plane-parallel slab is thus 
dz = H ds, with \x = cos 6 and where 8 is the angle between z 
and s (see Fig. lA.U . Equation lA.ll is thus written as 



fl— = S v - Iv 

dT 



where the differential optical depth is given by gas and dust 
contributions, dr = a v dz, and S v = j v /a v is referred to as the 
source function. Continuum emissivity from dust is assumed to 
be thermal and given by 



•dust dust d / T" \ 

j v = a v B v (T d ) 



(A.3) 



where B v is the Planck function at a given dust temperature, Tj, 
and dy* is computed from any of the dust mass absorption 
cross-sections available in the literature (e.g. Draine & Lee, 
1984, Ossenkopf and Henning, 1994). For practical purposes, 
the dust absorption coefficient is assumed to be constant in all 
the passband around each line frequency. Hence, the number 
of dust continuum photons emitted per second in a given cell 
of material is (An /he) j v ius ' V m Av where V m is the cell volume 
and Av the considered passband in velocity units. Although the 
inclusion of dust almost does not affect molecular excitation in 
our work, it is included for consistency and for making predic- 
tions of higher frequency lines where it has larger effects. 

Molecular lines occur at discrete frequencies, vy, where i 
and j refer to upper and lower energy levels with «, and n ; 
[cm -3 ] populations respectively. Gas emission and absorption 
coefficients, as a function of velocity, are defined as 



he 



he 

= —(njB ji -n i B i j)cf> (A.4) 



where Bp, B/j, and Ay are the transition probabilities, or 
Einstein coefficients, for absorption, and for induced and spon- 
taneous emission respectively. We have assumed the same line 
Doppler profile (in velocity units) for emission and absorption 



1 



b y/n 



exp 



/ V + Vf ■ s \ 2 

I b / 



(A.5) 



and thus considered Gaussian Doppler microturbulent and ther- 
mal broadening characterized by the broadening parameter 
b 2 = v 2 urb + \ 2 h . Note that any arbitrary velocity field Vf can 
be included. Here we take the possibility of having a velocity 
field normal to the slabs, i.e. Vf=v/ (z). 

Generally speaking, the relative level populations of a con- 
sidered molecule m are determined by collisions with other 
molecules, and/or by radiative effects caused by the cosmic 
background and/or by the dust continuum emission. The partic- 
ular physical conditions, type of molecule and spectral domain 
will determine the dominant processes through the steady state 
statistical equilibrium equations 



j*' 



Cnl-X"^"-^! 



rot levels 



X »J ( A - 6 ) 



(A.2) 



where Cy and Rjj [s l ] are the collisional and radiative tran- 
sition rates between i and j levels. For the collisional rates of 
species m (CS, C 34 S, C I8 and HCS + ) we have considered 

C U = (H 2 ) n(H 2 ) + yPj (He) n(He) + y£ (H) n(H) (A.7) 

where [cm 3 s~'] are the temperature dependent collisional 
de-excitation rate coefficients of m with collisional partners H2 
or He. If unknown, excitation rate coefficients are computed 
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through detailed balance. For consistency with the PDR mod- 
eling we have estimated the collisional rates with H atoms (sim- 
ply by scaling from the He rates), since H may be the dominant 
partner in the diffuse regions. Radiative rates are 

Rij = A-ij + />';;./,, ; Rji = B j, J ji (A.8) 

where /y is the intensity of the radiation field integrated over 
solid angles and over the line profile. External illumination by 
cosmic background, dust continuum emission and line photons 
from different spatial regions contribute to 7y. Hence, the solu- 
tion of molecular excitation implicitly depends on the nonlocal 
radiation field, which obviously depends on level populations 
in many cloud points, Jy is explicitly computed in the Monte 
Carlo approach, and thus, the RT-excitation problem is solved 
iteratively until desired convergence in some physical parame- 
ter (generally the level populations) is achieved. LTE level pop- 
ulations at a constant fictitious T ex were used for the first iter- 
ation. In the case of CS modeling, T ror from the observational 
rotational-diagrams (Fig0 was used. 

The RT problem is then simulated by the emission of a de- 
termined number of model photons (both sides external illumi- 
nation, continuum and line photons) in a similar way to that 
originally described by Bernes (1979). Model photons repre- 
sent a large quantity of real photons randomly distributed over 
the line profile Ia31 and emitted at random cloud positions and 
directions. Each model photon is followed through the differ- 
ent slabs until it escapes the cloud or until the number of rep- 
resented real photons become insignificant. Note that the angle 
8 between the photon direction and the normal to the slabs re- 
mains constant in all the photon path. In spherical geometry the 
8 angle between the photon direction and the radial direction 
changes in each photon step and thus has to be computed re- 
peatedly. In addition, model photons sent in the cos 8^0 direc- 
tion in semi-infinite plane-parallel geometry will almost never 
escape the cloud. Thus, a minimum number of represented real 
photons is defined otherwise the photon is not followed any- 
more. In this way, the Monte Carlo simulation explicitly com- 
putes the induced emissions caused by the different types of 
model photons in all the slabs. At the end of the simulation, an 
averaged value for the By/y that independently accounts for 
external illumination, continuum emission and line emission is 
stored for every slab (2 S ij, m in Bernes formalism). Populations 
are then obtained in each slab by solving: 

»<■ + Z s ^ + Q j ] = Z n ' [ f. Z Si j-»< + c fi ] (A - 9) 

A reference field for all types of model photons was included 
to reduce the inherent random fluctuations (i.e. the variance) 
of any Monte Carlo simulation (see Bernes 1979, for details). 
When a prescribed convergence in level the populations is 
reached, the total (gas+dust) source function is completely de- 
termined and the emergent intensity can be easily computed by 
integrating Ea. IA.2l 

For spherical geometry, the code was successfully bench- 
marked against two test problems, the Bernes' CO cloud 
(Bernes 1979) and the HCO + collapsing cloud, problem 2a 
of the 1999 Leiden benchmark (van Zadelhoff et al. 2002). 



In the case of plane-parallel geometry, several thermalization 
tests for the CS excitation (without dust emission) were suc- 
cessfully performed. Excitation temperatures as a function of 
the normal coordinate to slabs z are shown for CS 7=1-0, 2- 
1, 3-2 and 5-4- transitions in Fig. IA.2I Model parameters are 
Tfc=20 K, n(H 2 )=10 5 cm" 3 and^(CS)=7xlO- 9 . Different exci- 
tation conditions are considered. Upper model: Ay/Cy = and 
T ex is correctly thermalized to T*,-„ (LTE). Middle model: colli- 
sional rates from Lique et al., (filled squares; 2006) and Turner 
et al. (empty circles; 1992) and resulting non-LTE excitation 
(see section l5Tl . As noted by Lique et al., their new collisional 
rates produce larger excitation temperatures, especially as / 
increases. For the Horsehead physical conditions this implies 
that the estimated densities and/or abundances with Turner et 
al. collisional rates are ~10% larger. Lower model: collisional 
excitation neglected and T ex is radiatively thermalized to the 
cosmic background temperature at Tj s =2.7 K. 

A.2. A model for an edge-on cloud with inclination 

In order to benchmark the spatial distribution of the PDR code 
abundance predictions with our interferometric line observa- 
tions, we now can use the simple model described above to 
compute the synthetic spectrum of a required molecule. To do 
that, the PDR code output was used as an input for the RT 
calculation. In particular, the density profile (both n(H 2 ) and 
n(H)), temperature profile (both T^ and T^) and species abun- 
dance are carefully interpolated from the PDR spatial grid out- 
put. In practice, the slab discretization for the RT calculation 
has to be precise enough to sample the abundance, density and 
temperature variations at the edge of the PDR (where most of 
the changes occur). In most RT computations, ~50 slabs were 
judged to give satisfactory sampling of the PDR variations. 
For an edge-on configuration, after a Monte Carlo simulation, 
RT equation lA. 21 was integrated in a grid of different lines of 
sight (similar to impact parameters in spherical geometry) as 




Fig. A.3. Adopted geometry for a plane-parallel PDR inclined 
by an small angle <p relative to the line of sight s. In this sketch, 
z denotes the normal direction to the slabs and also the UV 
illumination direction. 
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Fig. A.2. Thermalization tests for a plane-parallel cloud illuminated by the cosmic background at both surfaces (z/z m ax=0,l). 



depicted in Fig. IA.3l Lines of sight can be inclined by an angle 
<p respect to the slabs normal (ds = dz/sin <p). Therefore, the 
maximum integration path is Idepthlcos tp where Idepth is the as- 
sumed cloud spatial depth. To produce a synthetic map, results 
are then convolved in a grid of cloud points with an angular 
resolution characterized by a gaussian with hpbw equal to that 
of the synthesized interferometric beam. 

References 

Abergel, A., Bernard, J. P., Boulanger, F. et al. 2002, A&A, 389, 239. 
Abergel, A. et al. 2003, A&A, 410, 577. 
Anthony-Twarog, B. J. 1982, A&J, 87, 1213. 

Asplund, M., Grevesse, N., & Sauval, A. J. 2005 in ASP Conf. Ser. 
336, Cosmic Abundances as Records of Stellar Evolution and 
Nucleosynthesis, ed. F. N. Bash & T. G. Barnes, 25 

Batrla, W. & Wilson, T. L. 2003, A&A, 408, 231 

Bernes, C. 1979, A&A, 73, 67. 

Bogey, M., Demuynck, C, & Destombes, J. L. 1981, 

Chem.Phys.Letters, 81, 256 
Boogert, A. C. A. et al. 2000, A&A, 360, 683. 
Cernicharo, J. & Guelin, M. 1987, A&A, 176, 299. 
Dartois, E. 2005, Space Science Reviews, 119, Issue 1-4, p. 293-310 
Draine, B. T. 1978, ApJS, 36, 595. 
Draine, B. T. & Lee, H.M. 1984, ApJ, 285, 89. 
Emerson, D. T, & Graeve, R. 1988, A&A, 190, 353. 
Falgarone, E., Verstraete, L., Pineau Des Forets, G. & Hily-Blant, P. 

2005, A&A, 433, 997. 

Falgarone, E., Pineau Des Forets, G., Hily-Blant, P. & Schilke, P. 

2006, A&A, in press. 

Flower, D. R., & Pineau des Forets, G. 1998, MNRAS, 297, 1 182. 
Flower, D. R. 2001, JPhB, 34, 2731 

Frerking, M. A., Wilson, R. W., Linke, R. A., & Wannier, P. G. 1980, 
ApJ, 240, 65. 

Garcfa-Rojas, J., Esteban, C, Peimbert, M., Costado, M. T., Rodrguez, 

M., Peimbert, A. & Ruiz, M. T. 2006, MNRAS, 368, 253 
Gerin, M. et al. 1997, A&A, 318, 579. 

Gibb, E. L., Whittet, D. C. B., Boogert, A. C. A., & Tielens, A.G.G.M. 

2004, ApJS, 151,35. 
Goicoechea, J.R. 2003, Ph. D. Thesis, Universidad Autonoma de 

Madrid, September 2003. 
Goldsmith, P.F., & Langer, W.D. 1999, ApJ, 517, 209. 
Gonzalez-Alfonso, E. & Cernicharo, J. 1993, A&A, 279, 506. 
Graedel, T. E., Langer, W. D., & Frerking, M. A. 1982, ApJS, 48, 321. 
Graf, U. U„ Genzel, R., Harris, A. I., Hills, R. E., Russell, A. P. G. & 

Stutzki, J. 1990, ApJ, 358, L49. 
Gueth, F, Guilloteau, S. & Bachiller, R. 1996, A&A, 307, 891-897. 



Habart, E., Abergel, A., Walmsley, C. M., Teyssier, D. & Pety, J. 2005, 

A&A, 437, 177-188. 
Hily-Blant, P., Teyssier, D., Philipp, S. & Gusten, R. 2005, A&A, 440, 

909. 

Hosokawa, T. & Inutsuka S. 2005a, ApJ, 623, 917. 
Hosokawa, T. & Inutsuka S. 2005b, astro-ph/0511165 
Howk, J. C, Sembach, K. R., & Savage, B. D. 2006 2006, ApJ, 637, 
333. 

Irvine, W. M., Schloerb, F. P., Hjalmarson, A. & Herbst, E. 1985, 

Protostars and planets II (A86-12626 03-90). Tucson, AZ, 

University of Arizona Press, 1985, 579-620 
Jansen, D. J., Spaans, M., Hogerheijde, M. R., & van Dishoeck, E. F. 

1995, A&A, 303, 541. 
Joulain, K., Falgarone, E., Pineau des Forets, G. & Flower, D. 1998, 

A&A, 340, 241. 

Le Bourlot, J., Pineau Des Forets, G, Roueff, E., & Flower, D. R. 

1993, A&A, 267, 233. 

Le Petit, F, Nehme, C, Le Bourlot, J. & Roueff, E. 2005, ApJS, 

preprint, doi:10.10861'503252\ 
Lique, F, Spielfiedel, A. & Cernicharo, J. 2006, A&A, 451, 1125 
Lis, D. C, & Schilke, P. 2003, ApJ, 597, L145. 
Lucas, R., & Liszt, H. 1998, A&A, 337, 246. 
Lucas, R., & Liszt, H. 2002, A&A, 384, 1054. 
Martin-Hernandez et al. 2002 2002, A&A, 381, 606 
Martin, S., Martin-Pintado, J., Mauersberger, R., Henkel, C, & 

Garcfa-Burillo, S. 2005, ApJ, 620, 210 
Millar, T. J. & Herbst, E. 1990, A&A, 231, 466. 
Montaigne, H., Geppert, W. D., Semaniak, J., et al. 2005, ApJ, 631, 

653. 

Monteiro, T. 1984, MNRAS, 210, 1 

Neufeld, DA., Wolfire, M.G., & Schilke, P. 2005, ApJ, 628, 60. 

Ohishi, M. & Kaifu, N. 1998, Chemistry and Physics of Molecules 
and Grains in Space. Faraday Discussions No. 109. The Faraday 
Division of the Royal Society of Chemistry, London, 1998, 205 

Ossenkopf, V. & Henning, Th. 1994, A&A, 291, 943. 

Palumbo, M.E., Geballe, T.R., Tielens, A.G.G.M. 1997, ApJ, 479, 
839. 

Pankonin, V., & Walmsley, CM. 1978, A&A, 64, 333. 
1973, ARA&A, 11, 51. 

Pety, J., Teyssier, D., Fosse, D., Gerin, M., Roueff, E., Abergel, A., 
Habart, E. & Cernicharo, J. 2005, A&A, 435, 885-899. 

Pety, J. 2005, in SF2A-2005: Semaine de l'Astrophysique Francaise, 
721 

Pound, M. W., Reipurth, B. & Bally, J. 2003, A&J, 125, 2108-2122. 

Reipurth, B., & Bouchet, P. 1984, A&A, 137, 1. 

Ruffle, D. P., Hartquist, T. W., Caselli, P. & Williams, D. A. 1999, 

MNRAS, 306, 691. 
Tieftrunk, A., Pineau des Forets, G., Schilke, P. & Walmsley, C. M. 

1994, A&A, 289, 579. 



22 



J.R. Goicoechea et al.: Low sulfur depletion in the Horsehead PDR 



Teyssier, D., Fosse, D., Gerin, M., Pety, J., Abergel, A. & Roueff, E. 

2004, A&A, 417, 135-149. 
Turner, B. E., Chan, K.W, Green, S., & Lubowich, D. A. 1992, ApJ, 

399, 114. 

van der Tak, F.F.S., Boonman, A.M.S., Braakman, R., & van 

Dishoeck, E.F. 2003, A&A, 412, 133 
van Dishoeck, E. F. 1988, Rate Coefficients in Astrochemistry. 

Editors, T.J. Millar, DA. Williams; Publisher, Kluwer Academic 

Publishers, Dordrecht, Boston, 49. 
van Zadelhoff, G.-J., Dullemond, C. P., van der Tak, EES. et al. 2002, 

A&A, 395, 373. 

Wakelam, V., Caselli, P., Ceccarelli, C, Herbst, E. & Castets, A. 2003, 

A&A, 422, 159 
Wannier, P. G. 1980, ARA&A, 18, 399. 
Zhou, S„ Jaffe, D. T, Howe, J. E. et al. 1993, ApJ, 419, 190. 



